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We study vibrational energy transport in a quasi 1-D harmonic chain with both longitudi- 
nal and transverse vibrations. We demonstrate via both numerical simulation and theoretic 
analysis that for 1-D atomic chain connected by 3D harmonic springs, the coefficient of heat 
conduction changes continuously with its lattice constant, indicating the qualitative differ- 
ence from the corresponding 1-D case where the coefficient is independent of the lattice 
constant. 
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1. Introduction 

Heat conduction in low dimensional systems (less than three dimension) has been attract- 
ing increasing attention in the past decade. 1 From the fundamental point of view, one would 
like to know whether the fundamental transport theory for bulk material, such as the Fourier 
law of heat conduction, is still valid for such low dimensional systems. In fact, the question is 
not trivial, as there is still no rigorous proof available so far. From application point of view, 
it is an indispensable question to understand the heat conduction properties of the nanoscale 
materials before they are put into application. 2 

In order to understand the underlying physical mechanism of heat conduction in one 
dimensional (ID) systems, different lattice models with and without on-site (pinning) potential 
have been used, such as the Fermi-Pasta-Ulam model 3 without on-site potential, the Frenkel 
Kontorova model, 4 and the </> model with on-site potential. 5 With these models, the roles of 
anharmonicity and the on-site potential in heat conduction in ID systems have been nicely 
demonstrated. For example, in the system without on-site potential, a size-dependent thermal 
conductivity has been observed 3 due to a superdiffusive motion of phonons. 6 In the system 
with on-site potential, the scattering from on-site potential makes phonon transport diffusively, 
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leading to a size-independent thermal conductivity. Considering the fact that in the ID models, 
the lattice is restricted to longitudinal vibration only, i.e., the transverse motions have been 
completely ignored, we conclude that these ID models are too simple to be used for modelling 
heat conduction in realistic nanostructures, such as nanotube and nanowires. 

Indeed, the heat conduction behavior changes when the transverse vibrations are consid- 
ered. In a recent work, Wang and Li 7 considered a model of quasi ID chain connected by 
2D spring with a bending angle interaction. The lattices are allowed to vibrate in longitudi- 
nal as well as in one transverse direction. The results show that for a fixed lattice constant, 
the anomalous heat conduction coefficient n(N) changes from a logarithmic divergence with 
system size N, i.e, logN for large transverse coupling, to a 1/3 power law divergence, iV 1 / 3 , 
at intermediate coupling, and then to a 2/5 power law divergence, N 2 / 5 at low temperatures 
and weak coupling. The results are mainly due to the mode coupling between the longitudinal 
modes and transverse modes. 

In this paper, we study the quasi ID chain with motion in both the longitudinal direction 
and the two transverse directions. Namely, the lattices can vibrate in all three directions, 
which is more closer to real nano-scale quasi ID systems. For simplicity, we call our system 
a 3D harmonic chain. To be more practical, we will not going to discuss the thermodynamic 
limit (length goes to infinity) as most of the nanoscale systems are of finite length. Instead, 
we shall focus on the effect of the lattice constant. The model discussed by Wang and Li 7 can 
be considered as a simplified polymer chain, while the model studied in our current paper can 
be considered as simplified nanotube or nanowire model. 

We shall demonstrate soon that, in ID case, the lattice constant does not play any role 
when the atoms can be considered as moving around its equilibrium; while in 3D case, the 
situation is totally different and the lattice constant will influence the heat conduction. It goes 
back to the ID case when the lattice constant a — > oo. 

2. Model and numerical results 

The 3D harmonic chain is coupled with the Nose-Hoover thermostats 8 on the first and 
last particle, keeping them at temperature Th and 7], respectively. The Hamiltonian is 



where a is the lattice constant, r^ = (ia + Xi,yi, Zi), and Xi, yi, Z{ are the displacements from 
the equilibrium position (ia, 0,0). The motion of the particles satisfy the canonical equations 
Qiw = dH/dpi w and pi w = —dH/dqi w with w = x,y, z, respectively, and i = 2, 3, ■ ■ ■ , N — 1. 
The left and right heat baths are described respectively by the equations, 




(1) 



ih = ^(<?L + Qly + Qlz) ~ 1, 
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il = ^f (.$fx + <lNy +q 2 Nz) ~ 1. ( 2 ) 
The equations of the motion for the first and last particle are q\ w = dH/dpi w , p\ w = 
-dH/dq lw - ChPiw, QNw = dH/dp Nw , and p Nw = -dH/dq Nw - iiPNw with w = x,y,z, 
respectively. Substituting Eq. (1) into the canonical equations we can get the expressions of 
Pix, Piy and pi z as follows 

Pix = (1 - -r— ){a + Ax i+ i) - (1 - ^-){a + Axi), 

«i+l &i 

Piy = {l-^-)Ay l+1 -{l-j)Ay u (3) 

Piz = (1 - -^—)Az i+1 - (1 - ^-)Az u 
di+i di 

where d i+1 = y/(a + Ax i+1 ) 2 + (Ay i+1 ) 2 + (Az 1+ i) 2 , Ax i+1 = x i+1 - x*, Ay i+1 = y i+1 - y u 
and Az i+ i = z i+1 - z { . 

Re-scaling q iw and p iw as q' iw = qi w ja and p' iw = Pi w /a with w = x,y, z, respectively, Eq. 
(3) becomes 

(1- J-)(l + A^ +1 )-(l-^)(l + A^), 



P'iy = (1 " ^-)A^ +1 - (1 - ^)Ay' h (4) 
a i+i 



p' lz = a - -J-)^+i - a - i)A^ 

where = y^l + Ax- +1 ) 2 + (Ay 4 ' +1 ) 2 + (Az' i+1 ) 2 . The heat bath Eqs. (2) is re-scaled as: 



& = gjv^i* + ^ + - 1, 

6 = ypWNx + Q'n V + Q'Nz) - !■ 



1 , -.2 -.2 -,2 

if 

1 / -.2 -,2 -,2 

where T)£ 2 = Thj/a 2 . This means that the equations pi^ and pN w for the boundary particles 
are invariant providing that the temperature of the two heat baths are re-scaled as mentioned 
above. This is interesting. It indicates that in our system, the model with twice smaller lattice 
constant is identical with the one with four times smaller temperature. 

For the case of fixed temperatures Th and I), the re-scaling of q\ w and q^w will cause 
different £/, and namely, the effect of lattice constant will show up through the heat baths 
and result in different behaviors of heat conduction. This can be understood like this. For 
a given Th and T}, the average of the amplitude of <ji w are fixed for the heat bath particles. 
Therefore, larger a means a relatively smaller q' iw . Making a Taylor expansion of d' i+1 with 
respect to Ax^, Ay-, and Az[, from Eq. (4) we see that p' ix is proportional to the first order 
of Ax^ and p' iy and p' iz are proportional to the second order of Ax^. Therefore, smaller q' iw or 
larger a will make the transverse motion negligible, indicating that we may observe a harmonic 
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motion for the larger a. Generally, with the decrease of a, the transverse motion will become 
more and more important and thus influence the heat conduction more. When q' iw reaches 
the magnitude of unity or of lattice constant, the heat conduction will be saturated, thus we 
may observe a platform for small enough a. 

3. Heat conductivity versus system size and lattice constant 

The temperature of the i'th lattice is defined as T(i) = ((pf x ) + (pf y ) + (p 2 z ))/3 and the heat 
flux along the chain in stationary state is defined as J = (pix^^ 1 + Piy ^fjj^ 1 + Piz^^r 1 ) 
where Vi+i^i = ^(1^+1 — r»| — a) 2 . From the Fourier law J = — kVT, the heat conductivity 
can be expressed as 



For fixed and 7), JN can be used to represent the k. For a momentum conservative ID 
lattice, it is found that the heat conduction is anomalous, 1 ' 3 ' 6 namely, the heat conductivity 
depends on the system size as 



where (3 varies from system to system. In this paper, we study the dependence of (3 on the 
lattice constant a. For comparison, we study two concrete cases: a relatively high temperature 
regime and a relatively low temperature regime. For the low temperature case, we take = 
0.07 and 7] = 0.05, thus the system temperature defined by T = (T h + 7])/2 is T = 0.06. 
For high temperature case, we take 7X = 0.3 and 7} = 0.2, thus the system temperature is 
T = 0.25 which is about four times higher than the previous one. The results are shown in 
Fig. 1(a) and (b). Our numerical simulations reveal that NJ increases with N by different 
scaling for different a and different temperatures T h / t . 

As we have pointed out early that, in the dimensionless Hamiltonian, the system for 
different a and different heat bath temperatures can be equivalent, therefore, in Fig 1(c) we 
draw (3 versus the quantity a/VT. 

It is quite obvious that two sets of data from Fig. 1(a) and (b) almost overlap with each 
other in Fig. 1(c). We have also checked other pairs of 7X and 7} ( such as 7\ = 0.1 and 
7] = 0.005) and found that they also overlap with each other. From Fig. 1(c) it is easy to 
see that (3 goes to unity for a/VT > 300, indicating the large lattice constant or very low 
temperature makes the 3D chain equivalent to a ID harmonic chain. Moreover, (3 drops with 
the decrease of a or increase of temperature, and reaches a platform for a/VT < 15, confirming 
the theoretic analysis of the previous section. 

We have checked the behaviors of the particles and found that the particles oscillate 
around their equilibrium positions and can be considered as local when a is relatively large 
and their oscillation can not be considered as local when a is in the range of platform. For 
the relatively small a, it is possible for the particle to move to the area with \xi\ > a. Figure 2 



K(N) = J(N-T)/(T h -T l ). 



(6) 



K 



(N) ~ N 13 , 



(7) 
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Fig. 1. NJ versus N for different temperatures, (a) T h = 0.07 and T, = 0.05, thus T = (T h +T t )/2 = 
0.06. (b) T h = 0.3 and T t = 0.2, thus T = (T h + T t )/2 = 0.25. (c) versus a/VT. 

shows the snapshots of the deviation of particle's positions from their equilibriums. It can be 
seen that the positions at i = 17 — 23 satisfy \xi\ > a(= 2), and the similar situation occurs in 
Fig. 2(b) for % = 12. This kind of non-local behavior causes the platform. High environment 
temperature causes large deviation of \x^\ from a, hence results in a larger platform as is shown 
in Fig. 1(c). 

4. Anomalous energy diffusion 

In order to confirm further the influence of lattice constant on the heat conduction in the 
3D harmonic chain, we study the diffusion for different a. First, we let the two heat baths have 
the same temperature and let the chain run enough time to achieve the equilibrium. Then we 
give a pulse of energy, which is larger than the energy of the equilibrium state, at the middle 
particle of the chain and investigate how the pulse of energy spreads to the other part of the 
chain. Figure 3 shows the result for 7} = Th = 0.2, N = 201, and the results are obtained by 
averaging over 10 3 realizations where the initial pulse of five times the average momentum is 
added at the particle i = 101, i.e., at t = we change the velocity of the middle particle to 5 
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Fig. 2. Snapshots of Xi, the deviation from the equilibrium position, for a = 2. (a) Th = 0.07 and 
Ti = 0.05; (b) Th = 0.3 and Ti = 0.2. 

times of its previous value. Figure 3(a) and (a') denote the snapshot at time t = 1, (b) and 
(b') at t = 3, (c) and (c') at i = 10, (d) and (d') at t = 30, and (e) and (e') at t = 100, and 
the left panels represent the situation of a = 2 and the right panels represent the situation of 
a = 128. 

Obviously, the left panels of Figure 3 are different from the right panels, reflecting the 
influence of the lattice constant. Ref. 6 gives an approach to measure the energy diffusion by 
the formula, <J 2 {t) = ^ E ^E(x E t)-E )dx ^ ' w ^ ere E(x,t) is the energy distribution at time t, 
xq is the position of initial energy pulse at t = 0, Eq is the energy of the equilibrium state 
with temperature T. However, in the 3-D case, this approach does not work very well because 
the energy E(x,t) fluctuates around Eq and makes the diffusion contributions in the integral 
of this equation cancel each other. 

A correct way is to make the diffusion contributions always positive, i.e., make the numer- 
ator in this equation positive at each position x. Ref. 9 points out recently that, if the pulse is 
very small and is initially localized at i = io, the energy diffusion can be given by 

v 2 (t) = Y J (i-*o) 2 h)(^ t )i (8) 

i 

where (5(2) (i,t) = (5vi(t)) 2 in the ID case. Here we extend this approach to the 3D case and 
take 

<y (2) (M) = {^ x {t)f + {8 Piy {t)f + {5 Plz {t)f 
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Fig. 3. Snapshots of (E(i)) at different times for N — 201, T — 0.2, and the initial pulse plioia:(0) = 
5pioiz(0), pi 1012,(0) = 5pioi v (0), plioiz(O) = 5pioia(0) is added at the particle i = 101 at t = 
where the left panels represent the situation of a = 2, the right panels represent the situation of 
a = 128, and (a) and (a') denote the snapshot at time t = 1, (b) and (b') at t = 3, (c) and (c') at 
t = 10, (d) and (d') at t = 30, and (e) and (e') at t = 100. 

+(5x i (t)) 2 + (5y i (t)) 2 + (6z i (t)) 2 . (9) 

That is, (J 2 (t) measures the total energy variation caused by the diffusion of the pulse. In 
numerical simulations, after the system arrives at the equilibrium we record the time as 
t = and copy all the variables Xj and pi to a new set of variables xlj and plj except the 
middle one plioi = l-5pioi where 1.5 chosen so as to make the initial pulse small enough. 
Hence, we have two set of systems (xj,pj) and (xlj, plj with xlj = Xj, plj = pj for i ^ iq 
and xlj = Xj , pl i() = 1.5pj at t = 0. As time goes on, the corresponding variables will 
gradually become different because of the diffusion of the pulse. Then, we can obtain a nonzero 
5pix{t) = P^ix(t) —Pix(t), 5xi(t) = xli(t) — Xi(t), and so on. In this way, we can find the relation 
cr 2 (i) ~ t a through Eqs. (9) and (8). Figure 4(a) shows the result for three typical values of a 
at T = 0.05, which are obtained by averaging over 10 3 realizations. The line with "squares" 
denotes the case of a = 2, the line with "circles" the case of a = 16, and the line with 
"triangles" the case of a = 128. We have the similar situations for other a. The slopes of the 
lines in Fig. 4(a) give the exponent a. 
Fig. 4(b) show the case of T = 0.2. 

For the same reason as in Fig 1, we also draw here a versus a/yT in 4(c). The two sets 



7/10 



J. Phys. Soc. Jpn. 



Full Paper 



3.0 n 



2.5 



2.0- 



1.0- 



0.54 



(a) 



T=0.05 



3.0 



* « 2.5 

/ 

/ .' 2.0 



a=2 

a=16 

a=128 



1.5 
1.0 
0.5 



(b) 



T=0.2 



a=2 

a=16 

a=128 



0.6 0.9 1.2 1.5 0.6 0.9 1.2 1.5 



lg(t) 



lg(t) 



2 


0- 


1 


.9- 


1 


8- 


1 


7- 


s 1 


.6- 


1 


5- 


1 


4- 


1 


3- 



(c) 



.o 



-O- T=0.05 
-•- T=0.2 



10 100 



Fig. 4. Three typical situations of a versus t for T = 0.05 (a) and T = 0.2 (b). We have taken an 
ensemble average over 10 3 realizations, (c) a versus a/VT for T = 0.05 and T = 0.2. 

of date overlap quite well. More interestingly, the Fig. 4(c) clearly indicates that 

a = A + B log 10 (a/Vr) , (10) 

where A ~ 1.138 and -B ~ 0.325 are two constants. Why the exponent of superdiffusion 
is related to the lattice constant a and temperature T in such a relation deserves further 
investigation. 

5. Conclusions and discussions 

In summary, we have studied heat conduction in a quasi 1-D system with 3D vibration. 
For finite length (around 1,000 units or, for a typical crystal atom chain, 100 nanometers), 
the heat conduction is anomalous and energy diffusion is superdiffusive. Moreover, we have 
demonstrated that both the anomalous heat conduction exponent f3 and the anomalous diffu- 
sion exponent a depends on the dimensionless quantity, aj\[T . In particular, we found that 
a = A + B\n(a/VT) for reasons to be further investigated. 
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We should emphasize that the results obtained in this paper are for a finite size system, 
which is more appropriate as the nanoscale systems are always of finite size. Therefore, the 
obtained results may shed lights in understanding heat conduction in nanoscale structures, 
especially in networked structure. 10 As the quasi ID systems such as the nanowires and 
nanotubes can be fabricated easily these days, it is therefore of great interest to do the 
measurement on the size dependent thermal conductivity of these systems and to compare 
with the results presented. 
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